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GENERAL  FRAMEWORK  FOR  ACHIEVING  TEXTBOOK  MULTIGRID  EFFICIENCY: 

ONE-DIMENSIONAL  EULER  EXAMPLE 


JAMES  L.  THOMAS*,  BORIS  DISKIN+,  ACHI  BRANDT* ,  AND  JERRY  C.  SOUTH,  JR.§ 

Abstract.  A  general  multigrid  framework  is  discussed  for  obtaining  textbook  efficiency  to  solutions  of  the 
compressible  Euler  and  Navier-Stokes  equations  in  conservation  law  form.  The  general  methodology  relies  on 
a  distributed  relaxation  procedure  to  reduce  errors  in  regular  (smoothly  varying)  flow  regions;  separate  and 
distinct  treatments  for  each  of  the  factors  (elliptic  and/or  hyperbolic)  are  used  to  attain  optimal  reductions 
of  errors.  Near  boundaries  and  discontinuities  (shocks),  additional  local  relaxations  of  the  conservative 
equations  are  necessary.  Example  calculations  are  made  for  the  quasi-one-dimensional  Euler  equations;  the 
calculations  illustrate  the  general  procedure. 

Key  words,  textbook  multigrid  efficiency,  distributed  relaxation,  Euler  equations 

Subject  classification.  Applied  and  Numerical  Mathematics 

1.  Introduction.  Computational  fluid  dynamics  (CFD)  has  become  an  integral  part  of  the  aircraft 
design  cycle  because  of  the  availability  of  faster  computers  with  more  memory  and  improved  numerical 
algorithms  and  physical  models.  More  impact  is  possible  if  reliable  methods  can  be  devised  for  off-design 
performance,  generally  associated  with  unsteady,  separated,  vortical  flows  with  strong  shock  waves.  Such 
computations  demand  significantly  more  computing  resources  than  are  currently  available. 

The  current  Reynolds- averaged  Navier-Stokes  (RANS)  solvers  with  multigrid  algorithms  require  on  the 
order  of  1500  residual  evaluations  to  converge  the  lift  and  drag  to  one  percent  of  their  final  values,  even  for 
relatively  simple  wing-body  geometries  at  transonic  cruise  conditions.  It  is  well  known  for  elliptic  problems 
that  solutions  can  be  attained  optimally  by  using  a  full  multigrid  (FMG)  process  in  far  fewer  (on  the  order  of 
2  to  4)  residual  evaluations.  A  multigrid  method  is  defined  by  Brandt  [2,  3,  4]  as  having  textbook  multigrid 
efficiency  (TME)  if  the  solutions  to  the  governing  system  of  equations  are  attained  in  a  computational  work 
that  is  a  small  (less  than  10)  multiple  of  the  operation  count  in  the  discretized  system  of  equations.  Thus, 
operation  count  may  be  reduced  by  several  orders  of  magnitude  if  TME  can  be  attained  for  the  RANS 
equations. 

State-of-the-art  multigrid  methodologies  for  large-scale  compressible  flow  applications  use  a  block-matrix 
relaxation  and/or  a  pseudo-time- dependent  approach  to  solve  the  equations;  significant  improvements  have 
been  demonstrated  with  multigrid  approaches,  but  the  methods  are  not  optimally  convergent.  The  RANS 
equation  sets  are  systems  of  coupled  nonlinear  equations  which  are  not,  even  for  subsonic  Mach  numbers, 
fully  elliptic,  but  contain  hyperbolic  partitions.  The  distributed  relaxation  approach  of  Brandt  [2,  3]  de¬ 
composes  the  system  of  equations  into  separable,  usually  scalar,  factors  that  can  be  treated  with  optimal 
methods.  Several  years  ago,  an  investigation  was  started  to  extend  this  approach  to  large-scale  applications; 

*  Computational  Modeling  and  Simulation  Brandi,  Mail  Stop  128,  NASA  Langley  Research  Center,  Hampton,  Virginia  23681 
(email:  j  .  1 .  thomasQlarc .  nasa .  go v) . 

^Institute  for  Computer  Applications  in  Science  and  Engineering,  Mail  Stop  132C,  NASA  Langley  Research  Center,  Hampton, 
Virginia  23681  (email:  bdiskinQicase .  edu).  This  research  was  supported  by  the  National  Aeronautics  and  Space  Administration 
under  NASA  Contract  No.  NAS  1-97046  while  the  author  was  in  residence  at  the  Institute  for  Computer  Applications  in  Science 
and  Engineering  (ICASE),  NASA  Langley  Research  Center,  Hampton,  Virginia  23681-2199. 

*The  Weizmann  Institute  of  Science,  Rehovot  76100,  Israel  (email:  achiQwisdora.weizmann.ac.il). 

§  Williamsburg,  Virginia  23185 


1 


at  that  time,  several  TME  demonstrations  for  incompressible  simulations  had  been  completed  and  Ta’asan 
had  shown  promising  results  for  the  subsonic  Euler  equations  [18].  Progress  has  been  shown  in  extending 
the  methodology  to  viscous  compressible  flow  applications  [19]  and  to  compressible  Euler  equations  using 
a  compact  differencing  scheme  [17].  Further  incompressible  flow  applications  have  been  made,  including 
complex  geometries  [14]  and  high-Reynolds-number  viscous  flow  in  two  [20]  and  three  dimensions  [13].  The 
progress  and  remaining  barriers  in  TME  for  the  equations  of  fluid  dynamics  were  summarized  in  [4]. 

The  purpose  of  this  paper  is  to  discuss  the  general  framework  expected  to  be  required  for  large-scale 
compressible  flow  applications.  The  quasi-one-dimensional  Euler  equations  are  solved  to  illustrate  the  frame¬ 
work.  Fully  subsonic  and  supersonic  applications,  as  well  as  transonic  applications  with  a  captured  shock, 
are  shown. 

2.  General  Framework.  The  viscous  compressible  equations  for  the  time-dependent  conservation  of 
mass,  momentum,  and  energy  can  be  written  as 


ftQ  +  R  =  0,  (2.1) 

where  the  conserved  variables  are  Q  =  (pypu1  pv,  pw,pE)T,  representing  the  density,  momentum  vector, 
and  total  energy  per  a  unit  volume,  and  R(Q)  is  the  spatial  divergence  of  a  vector  function  representing 
convection  and  viscous  and  heat  transfer  effects.  In  general,  the  simplest  form  of  the  differential  equations 
corresponds  to  nonconservative  equations  expressed  in  primitive  variables,  here  taken  as  the  set  composed  of 
velocity,  pressure,  and  internal  energy,  q  =  (u,  v,w,p,  e)T.  These  equations  are  found  readily  by  transforming 
the  time-dependent  conservation  equations  to  time-dependent  primitive  variable  equations.  Similarly,  a  set  of 
nonconservative  correction  equations  can  be  derived,  with  a  right  hand  side  vector  composed  of  a  combination 
of  the  conserved  residual  terms,  given  as 


L«5q  =  -^R.  (2.2) 

In  Eq.  (2.2),  is  the  Jacobian  matrix  of  the  transformation  and  the  correction  Sq  =  qn+1  -  qn,  where  n 
is  an  iteration  counter.  For  steady-state  equations,  the  time  derivative  is  dropped.  At  the  discrete  level,  the 
right  side  of  the  correction  equation  (2.2)  is  a  combination  of  conservative-discretization  residuals  while  the 
left  side  is  the  principal  linearization  of  the  nonconservative  operator. 

Note  this  is  not  a  Newton  linearization;  only  the  principal  terms  in  a  linearization  of  R  are  retained. 
The  principal  terms  are  those  terms  that  make  a  major  contribution  to  the  residual  per  a  unit  change  in 
q.  The  principal  terms  thus  generally  depend  on  the  scale,  or  mesh  size,  of  interest.  For  a  scalar  equation, 
the  discretized  highest  derivative  terms  are  principal  on  grids  with  small  enough  mesh  size  h.  In  deriving 
the  principal  linearization  for  high-Reynolds-number  simulations,  it  is  essential  to  consider  both  inviscid 
and  viscous  scales  —  the  inviscid  scales  dominate  over  most  of  the  flow  field  and  the  viscous  scales  are 
important  in  the  thin  viscous  layers  near  bodies  and  in  their  wakes.  Note  that,  for  a  discretized  system  of 
differential  equations,  such  as  R  =  0,  the  principal  terms  are  those  that  contribute  to  the  principal  terms  of 
the  determinant  of  the  matrix  operator  The  coefficients  of  the  principal  terms  in  L  are  evaluated  from 
the  current  approximation. 

The  principal  linearization  is  applied  to  the  correction  equation  based  on  a  nonconservative  approxi¬ 
mation.  Thus,  we  expect  the  correction  to  be  good  away  from  discontinuities  (shocks,  slip  lines)  in  the 
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flow  field.  It  is  in  these  regular  (smoothly  varying)  flow  regions  that  we  apply  distributed  relaxation.  The 
distributed  relaxation  method  replaces  by  M5w  so  that  the  resulting  matrix  L  M  becomes  a  diagonal  or 
lower  triangular  matrix,  as 


LM6w  =  (2-3) 

The  diagonal  elements  of  L  M  are  composed  ideally  of  the  separable  components  of  the  determinant  of  the 
matrix  L  and  represent  the  elliptic  or  hyperbolic  factors  of  the  equation.  In  [3,  8],  the  Sw  variables  were 
termed  as  “ghost  variables,”  because  they  need  not  explicitly  appear  in  the  calculations. 

The  distributed  relaxation  approach  yields  fast  convergence  for  both  steady  and  unsteady  simulations 
if  the  constituent  scalar  diagonal  operators  in  LM  are  solved  with  fast  methods.  The  approach  can  be 
applied  to  quite  general  equations;  a  set  of  matrices  M  has  been  derived  in  [3]  to  provide  a  convenient  lower 
triangular  form  for  the  compressible  and  incompressible  equations  of  fluid  dynamics  (including  a  variable 
equation  of  state). 

For  the  compressible  Euler  equations,  the  scalar  factors  constituting  the  main  diagonal  of  L  M  are  con¬ 
vection  and  full-potential  operators.  An  efficient  solver  for  the  former  can  be  based  on  downstream  marching, 
with  additional  special  procedures  for  recirculating  flows  [7,  8,  21];  the  latter  is  a  variable  type  operator,  and 
its  solution  requires  different  procedures  in  subsonic,  transonic,  and  supersonic  regions.  In  deep  subsonic 
regions,  the  full-potential  operator  is  uniformly  elliptic  and  therefore  standard  multigrid  methods  yield  opti¬ 
mal  efficiency.  When  the  Mach  number  approaches  unity,  the  operator  becomes  increasingly  anisotropic  and, 
because  some  smooth  characteristic  error  components  cannot  be  approximated  adequately  on  coarse  grids, 
classical  multigrid  methods  severely  degrade.  The  characteristic  components  are  those  components  that  are 
much  smoother  in  the  characteristic  directions  than  in  other  directions  [4,  11,  10].  In  the  deep  supersonic 
regions,  the  full-potential  operator  is  uniformly  hyperbolic  with  the  stream  direction  serving  as  the  time¬ 
like  direction.  In  this  region,  an  efficient  solver  can  be  obtained  with  a  downstream  marching  procedure. 
However,  this  procedure  becomes  problematic  for  the  Mach  number  dropping  towards  unity,  because  the 
Courant  number  associated  with  the  downstream  marching  procedure  is  large.  Thus,  a  special  procedure 
is  required  to  provide  an  efficient  solution  for  transonic  regions.  This  local  procedure  [5,  6,  9]  is  based  on 
piecewise  semicoarsening  and  some  rules  for  adding  dissipation  at  the  coarse-grid  levels. 

Boundaries  introduce  some  additional  complexity  in  distributed  relaxation.  The  determinant  of  LM  is 
usually  higher  order  than  the  determinant  of  L.  Thus,  as  a  set  of  new  variables,  Sw  would  generally  need 
additional  boundary  conditions.  In  relaxation,  because  the  ghost  variables  can  be  added  in  the  external 
part  of  the  domain,  it  is  usually  possible  to  determine  suitable  boundary  conditions  for  <5w  that  satisfy  the 
original  boundary  conditions  for  the  primitive  variables.  Examples  are  given  in  [20]  in  incompressible  flow 
for  entering  and  no-slip  boundaries.  However,  to  construct  such  a  remedy  may  be  difficult  and/or  time- 
consuming  in  general.  In  addition,  enforcing  these  boundary  conditions  causes  the  relaxation  equations  to 
be  coupled  near  the  boundaries,  not  decoupled  as  they  are  in  the  interior  of  the  domain. 

Thus,  near  boundaries  and  discontinuities,  the  general  approach  [3,  4]  is  to  relax  the  governing  equations 
directly  in  terms  of  primitive  variables.  Several  sweeps  of  robust  (but  possibly  slowly  converging)  relaxation, 
such  as  Newton-Kacmarcz  relaxation,  can  be  made  in  this  region.  The  additional  sweeps  will  not  affect 
the  overall  complexity  because  the  number  of  boundary  and/or  discontinuity  points  is  usually  negligible  in 
comparison  with  the  number  of  interior  points. 

This  general  framework  is  used  below  to  solve  the  quasi-one-dimensional  Euler  equations  in  fully  sub- 
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sonic,  fully  supersonic,  and  transonic  (with  and  without  shock)  flow  regimes.  The  regular  flow  regions 
are  relaxed  with  distributed  relaxation.  Boundaries  and  shocks  are  treated  by  applying  local  relaxation 
—  corresponding  here  to  updates  through  a  direct  solution  of  an  approximate-Newton  linearization  of  the 
conservative  equations. 

In  all  cases,  an  FMG-1  algorithm  is  used;  at  each  level,  the  equations  are  solved  with  an  FV(2,1)  full 
approximation  scheme  (FAS)  [2,  20]  multigrid  cycle.  Six  levels  are  used  in  the  cycle  wherever  possible.  The 
total  operation  count  in  the  FMG-1  algorithm  is  equivalent  to  about  40  residual  evaluations.  This  somewhat 
excessive  operation  count  is  typical  for  one-dimensional  problems,  where  the  number  of  coarse-grid  points  is 
large  relatively  to  multidimensional  cases.  In  three-dismensional  problems,  the  expected  operation  count  is 
about  8  residual  evaluations. 

3.  Quasi-One-Dimensional  Equations.  The  quasi-one-dimensional  equations  express  the  conserva¬ 
tion  of  mass,  momentum,  and  total  energy  as 

dt(  Q)  +  R  =  0,  (3.1) 


where 


Q  =  Q  a  =  (p,pu,pE)T  a, 

R  =  dx(Fa)  +  S, 

and  a  —  <r(x)  is  the  area  term.  The  flux  F  and  the  source  term  S  are  defined  as 


(3.2) 

(3.3) 


pu  ^ 

F  =  |  pu 2  +  p 

puE  +  up  } 


(3.4) 


S  =  - 


\  da 

,  p  Tx- 
\  0 


(3.5) 


The  pressure  p,  internal  (thermal)  energy  e,  and  sound  speed  c  are  related  through  the  equation  of  state  as 


P={  7  -  1  )pe. 
e  =  E-u2/2, 
<?  =  t p/p, 


(3.6) 

(3.7) 

(3.8) 


and  7  is  the  ratio  of  specific  heats. 

A  discrete  conservative  upwind-biased  differencing  approximation  to  R  can  be  defined  for  first-  or  second- 
order  accuracy  as 


-  (0,Pj,0)T[(rj+i  -ffj-j]- 


(3.9) 
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Here,  a  finite- volume  discretization  is  used,  where  subscripts  j  +  ~  and  j  —  |  denote  the  right  and  left 
interfaces,  respectively,  of  the  cell  centered  at  location  j  and  the  cell  spacing  in  the  x-direction  is  h .  The 
flux-differencing  splitting  of  Roe  [15]  is  used  to  construct  the  interface  flux  FJ+i;  pertinent  details  are 
described  in  appendix  A.  We  consider  only  first-  or  second-order  accuracy  and  do  not  differentiate  between 
average  and  pointwise  values  of  Q.  The  area  distribution  is  defined  as  a(x)  =  1  -  0.8x(l  —  x).  For  all  of 
the  results  presented  below,  we  overprescribe  the  boundary  values  from  the  exact  solution  of  the  differential 
problem  at  the  cell  centers  that  lie  outside  of  the  domain  0  <  x  <  1. 

4.  Relaxation  Schemes.  Several  schemes  are  considered  below  for  relaxation  of  the  steady-state  equa¬ 
tions 


R,  =  0.  (4.1) 

The  schemes  are  all  written  in  delta  form,  so  that  the  correction  to  the  solution  is  solved  at  each  iteration, 
denoted  as  SQ  =  Qn+1  -  Qn. 

4.1.  Conservative  Relaxation  Scheme.  The  exact  Newton  linearization  of  the  discrete  conserva¬ 
tive  equations  generally  leads  to  an  overly  complicated  linear  system  to  solve.  In  practice,  a  quasi-Newton 
(approximate)  linearization  is  used  to  reduce  the  algebraic  complexity  of  the  construction  and/or  the  band¬ 
width  of  the  resulting  linear  system.  For  the  purpose  of  constructing  a  conservative  relaxation  scheme,  an 
approximate  but  conservative  linearization  of  the  operator  (3.9)  can  be  written  as 


p-A++a+A-](Q).  (4.2) 

The  operators  S~  and  <5+  denote  backward  and  forward  first-order-accurate  difference  operators,  respectively, 
and  A+  and  A-  denote  the  positive  and  negative  eigenvalue  contributions  to  the  similarity  matrices  (see 
appendix  A  and  [12]  for  references). 

The  conservative  relaxation  update  <5Q  is  computed  from 


K-  A+  +  <5+  A-](*Q)  =  --R j.  (4.3) 

°3 

The  residuals  in  the  right  side  of  Eq.  (4.3)  are  computed  for  the  target  discretization  (4.1).  The  relaxation 
procedure  requires  a  tridiagonal  equation  to  be  directly  solved: 


(-A+_1,|A|j,A7+1)(JQ)  =  (4.4) 

Barth  [1]  analyzed  this  relaxation  with  fixed-point,  analysis  and  showed  the  relaxation  is  nearly  as  good 
as  a  full  Newton  iteration.  An  underrelaxation  parameter,  e.g.,  a  pseudo-time  step,  can  be  introduced  in 
the  relaxation  scheme  to  improve  robustness  and  stability.  Extensive  computations  made  during  the  course 
of  this  investigation  verified  that  this  relaxation  method  can  be  used  throughout  the  domain,  including 
shocks  and  boundaries.  In  the  general  framework  described  and  used  in  this  paper,  this  relaxation  is  applied 
only  locally  to  reduce  residuals  near  boundaries  and  shocks  and  to  solve  the  coarsest  grid  equation.  In 
multidimensional  case,  any  stable  procedure  converging  for  the  primitive-variable  conservative  equations 
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(4.1)  (e.g.,  pseudo-time  marching  relaxation)  can  be  applied  for  local  relaxation.  The  large  operation  count 
per  point  required  for  convergence  in  this  procedure  does  not  affect  the  overall  complexity  of  the  algorithm, 
because  the  number  of  grid  points  in  the  local-relaxation  regions  is  assumed  to  be  negligible  in  comparison 
with  the  total  number  of  grid  points. 

4.2.  Nonconservative  Relaxation  Scheme.  A  nonconservative  relaxation  scheme  is  derived  from  a 
conservative  scheme,  Eq.  (4.3),  with  the  assumption  that  the  Jacobian  matrices  are  locally  constant, 

[A  +*-  +  A-«5+](<5Q)  =  --Rj.  (4.5) 

aj 

This  formulation  is  expected  to  be  a  good  approximation  to  Eq.  (4.3),  except  near  discontinuities.  The 
corresponding  tridiagonal  equation  is 


(-At, \A\j,  A7)(«5Q)  =  -  Ar ..  (4.6) 

Subtracting  the  two  tridiagonal  equations,  (4.4)  and  (4.6)  yields 

(-A+ x  +  A+0,AJ+1  -  AJJOSQ)  =0,  (4.7) 

where  the  coefficients - A+_j  +  At  and  A“+1  —  AJ  —  are  0(h)  small  on  sufficiently  fine  meshes  in 

regular  flow  regions.  Near  shocks,  the  nonconservative  linearization  used  in  Eq.  (4.5)  does  not  satisfy  an 
order  property  and  the  nonconservative  relaxation  scheme  would  not  be  effective.  These  expectations  were 
confirmed  through  computations  with  both  the  conservative  and  nonconservative  relaxation  schemes,  where 
update  Eqs.  (4.4)  and  (4.6)  were  solved  precisely.  Barth  [1]  analyzed  a  similar  nonconservative  relaxation 
scheme  with  fixed-point  analysis  and  showed  poor  performance  of  the  nonconservative  scheme  for  a  flow  with 
a  shock. 

4.3.  Distributed  Relaxation.  Away  from  boundaries  and  shocks,  the  linear  update  scheme  Eq.  (4.5) 
is  relaxed  with  distributed  relaxation.  Transforming  the  corrections  in  Eq.  (4.5)  to  primitive  variables 
q  —  (u,p,e)T  gives 


[A +6~  +  A-,5+](<5q)  =  -(J^-Ar;,  (4.8) 

with  a  right-side  vector  composed  of  a  combination  of  the  conserved  residual  terms.  The  Jacobian  matrices 
are  related  through  a  similarity  transformation  to  the  conservative  Jacobians: 


A  = 


dq  dQ 
dQ  dq  ’ 


Alternately,  by  definition,  the  primitive  variable  correction  equation  is 


(4.9) 


L<5q  =  — r^, 


(4.10) 
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where 


=  R-. 

v  <9Q  h  a3  3 


(4.11) 


The  elements  of  the  matrix  operator  L  are  given  in  appendix  B. 

As  follows  from  Appendices  A  and  B,  the  determinant  of  L  for  (u  >  0)  in  deep  subsonic  or  supersonic 
regions  is 


det(L)  =  [( u 2  -  (?)  Sxx]uSx  , 


where 


(4.12) 


Sx  S +  for  fully  subsonic  flow, 
S~  S~  for  fully  supersonic  flow. 


(4.13) 


The  first  term  in  Eq.  (4.12)  represents  an  approximation  to  the  full-potential  operator  and  the  second 
term  represents  an  upwind  approximation  for  the  convection  of  entropy.  In  subsonic  flow,  Sxx  is  a  central 
discretization  associated  with  the  ellipticity  of  the  full-potential  operator.  In  fully  supersonic  flow,  Sxx  is 
a  one-sided  or  upwind-biased  approximation  in  accordance  with  the  hyperbolic  nature  of  the  full-potential 
operator.  For  first-order  upwind  differencing,  the  full-potential  factor  is  a  3-point  operator;  for  second-order 
upwind-biased  differencing,  the  factor  is  a  7-point  operator  (see  appendix  B). 

In  this  one-dimensional  case,  the  first  two  equations  in  the  matrix  operator  L  are  uncoupled  from  the 
third.  Thus  we  need  only  consider  distributed  relaxation  of  the  first  two  equations;  the  third  equation  can  be 
solved  for  St  in  its  primitive  variable  form  once  Su  and  Sp  are  found  from  distributed  relaxation.  Denoting 
the  upper  2x2  block  of  L  as  L,  an  obvious  choice  for  M  is  the  matrix  of  cofactors  of  L,  which  gives  a 
diagonal  matrix  for  LM  with  the  full-potential  factors  along  the  diagonal,  i.e., 


LM  = 


{u2-c?)Sxx 

0 


°  v 

{u2  -  c2)  6XX  ) 


(4.14) 


4.4.  Boundary  Treatment.  As  pointed  out  by  Sidilkover  [17],  there  is  a  connection  between  dis¬ 
tributed  relaxation  and  the  characteristic  variables.  For  subsonic  flow,  a  Jacobi  update  to  the  operators  in 
LM  is  equivalent  to  a  Kacmarcz  relaxation  of  the  characteristic  equations,  cast  in  terms  of  a  characteristic 
combination  of  corrections  as 


((u  4-  c)Sx  ( pcSu  +  Sp)  ]  _  I  Pcr i  +  r2 
(u  -  c)6+  (pcSu  -  Sp)  J  y  per  i  - 

With  first-order  differencing,  an  equal  and  opposite  correction  to  the  downstream  propagating  characteristic 
variable  pcSu  -I-  Sp  is  sent  in  Kacmarcz  relaxation  of  the  first  equation  to  the  local  cell  j  and  to  the  upstream 
cell  j  -  1;  the  correction  to  cell  j  is  exactly  one  half  of  that  found  from  a  point  Jacobi  relaxation  of 
the  corresponding  characteristic  equation.  Likewise,  an  equal  and  opposite  correction  to  the  upstream 
propagating  characteristic  variable  pcSu  —  Sp  is  sent  in  Kacmarcz  relaxation  of  the  second  equation  to  the 
local  cell  j  and  to  the  downstream  cell  j- hi.  With  a  local  update  of  the  solution  at  each  point,  the  smoothing 


(4.15) 
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Table  1 

The  discretization  errors  inp  at  convergence  and  the  relative  Li-norm  errors  after  the  FMG-1  algorithm  for  fully  subsonic 

flow. 


h 

IMI  :  P 

Ji£ii 

1-IMI  ■„ 
«*ll  • P 

1/32 

0.2654xl0-3 

<0.01 

1/64 

0.6341xl0-4 

<0.01 

1/128 

0.1547xl0~4 

<0.01 

1/256 

0.3817xl0-5 

<0.01 

rate  should  be  equivalent  to  Gauss-Seidel  relaxation  of  the  standard  3-point  Laplacian  operator.  However, 
the  relaxation  at  cells  adjacent  to  the  boundary  is  inconsistent  with  the  boundary  conditions  of  the  original 
problem;  for  example,  characteristic  boundary  conditions  would  prescribe  pcSu  4-  Sp  and  pcSu  -  Sp  as  zero 
at  the  upstream  and  downstream  boundaries,  respectively. 

This  example  demonstrates  the  general  result  mentioned  earlier:  that  distributed  relaxation  needs  to  be 
modified  near  the  boundaries.  A  simple  boundary  specification  for  the  ghost  variable  <5w  =  (Sw1,Sw2)T  can 
be  found  in  this  case;  enforcing  a  Neumann  condition  on  Swi  ( Sw2  )  and  a  Dirichlet  condition  on  Sw 2  (Swi 
)  at  inflow  (outflow)  in  subsonic  relaxation  is  consistent  with  characteristic  boundary  conditions  in  terms 
of  primitive  variables.  We  instead  use  the  more  general  formulation  and  apply  another  local  conservative 
relaxation  procedure  at  the  boundary  points  to  reduce  the  residuals.  Distributed  relaxation  is  then  applied 
in  the  interior  with  simple  Dirichlet  conditions  for  the  ghost  variables  Sw. 

4.5.  Relaxation  Sequence.  The  general  solution  procedure  was  to  reduce  the  local  residuals  at  least 
two  orders  of  magnitude  at  the  inflow  boundary,  the  outflow  boundary,  and  then  at  the  shock  region  —  both 
before  and  after  sweeping  the  entire  domain  with  distributed  relaxation.  The  local  conservative  relaxation  is 
made  over  the  first  two  cells  adjacent  to  the  boundary.  In  the  shock  region,  the  local  conservative  relaxation 
was  applied  to  nine  cells  centered  on  the  upstream  (supersonic)  side  of  the  shock.  A  pseudo-time  step 
corresponding  to  a  Courant  number  of  100  was  used  as  an  underrelaxation  factor  in  the  local  conservative 
relaxations;  otherwise  purely  steady-state  equations  were  relaxed.  In  distributed  relaxation,  the  variables 
(8wi,8w2)t  were  explicitly  used  in  the  implementation;  they  were  first  relaxed  over  the  interior  cells,  then 
distributions  to  (Ju,  8p)T  were  made,  and  then  the  corrections  to  Se  were  computed.  As  a  check,  single-grid 
computations  verified  that  the  convergence  of  the  error  per  relaxation  of  subsonic  equations  (4.1)  with  first- 
order  differencing  was  identical  to  that  found  with  a  Gauss-Seidel  relaxation  of  the  one-dimensional  3-point 
Laplacian  operator  with  Neumann  boundary  conditions  applied  at  one  end  of  the  domain. 

5.  Computational  Results.  Computational  results  are  shown  in  this  section  for  fully  subsonic,  fully 
supersonic,  and  transonic  flow  with  and  without  a  shock.  The  Mach  number  distribution  with  second-order 
accurate  discretization  on  a  grid  of  257  points  is  shown  in  Fig.  1  for  fully  subsonic,  fully  supersonic,  and 
transonic  flow  with  a  shock;  the  shock  location  is  specified  at  x  =  0.75.  There  is  no  limiting  applied  in 
the  transonic  cases  so  there  are  some  oscillations  at  the  shock.  In  this  one-dimensional  case,  the  solution 
can  easily  be  repaired  by  using  an  essentially  nonoscillatory  (ENO)  [16]  approach  (see  appendix  C)  but  the 
emphasis  of  this  work  is  on  the  solution  procedure  and  not  on  the  steady-state  results. 

5.1.  Fully  Subsonic  Flow.  The  smoothing  rate  of  distributed  relaxation  for  the  primitive  variable 
equations  is  expected  to  be  equal  to  the  worst  among  the  smoothing  rates  obtained  in  relaxation  of  the 
scalar  factors  of  the  determinant  of  L,  in  our  case  convection  and  the  full-potential  factor.  With  the  full- 
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Fig.  1.  Mach  number  distribution  computed  withoxtt  limiting  on  a  grid  of  257  points  for  fully  subsonic,  fully  supersonic, 
and  transonic  flow  with  a  shock  at  x  =  0.75. 


potential  operator  relaxed  (Gauss-Seidel)  and  the  convection  operator  solved  by  downstream  marching  from 
the  inflow  boundary,  the  convergence  rate  of  the  multigrid  cycle  per  relaxation  sweep  for  the  first-order 
accurate  discretization  should  be  close  to  the  smoothing  factor,  0.447,  of  Gauss-Seidel  relaxation  for  the 
one-dimensional  3-point  Laplacian  operator. 

Defect  correction  is  often  used  in  the  solution  of  higher  order  implicit  discretizations  to  reduce  the  arith¬ 
metic  operation  count  while  retaining  the  target  accuracy.  One  seemingly  possible  implementation  of  defect 
correction  is  to  use  first-order  accurate  discretizations  for  the  convection  operators  in  L.  The  corresponding 
distribution  matrix  M  consists  of  the  first-order  operators  as  well.  This  implementation,  however,  is  not  a 
good  idea  for  obtaining  good  smoothing  rates,  because  in  terms  of  high-frequency  components ,  a  low-order 
operator  and  a  corresponding  high-order  operator  do  not  necessarily  approximate  each  other  well.  Indeed, 
this  implementation  caused  a  slowdown  in  the  convergence,  and  the  smoothing  rate  corresponding  to  that 
of  the  Laplace  equation  was  not  attained. 

To  illustrate  this  phenomenon,  consider  the  Laplacian  operator  alone.  The  elliptic  operator  Sxx(wj) 
arising  from  differencing  the  primitive  variable  equations  with  the  second-order  upwind-biased  operator 
used  here  is  a  7-point  operator  (see  appendix  B).  Corrections  (Sw)j  to  the  current  approximate  solutions 
Wj  computed  in  Gauss-Seidel  defect-correction  relaxation  with  the  3-point  Laplacian  as  a  left-side  (driver) 
operator  can  be  found  from 


^2  [(Mj-i  -  2(<Mj]  =  (5’1) 

The  smoothing  factor,  computed  with  local  mode  analysis,  of  this  defect-correction  relaxation  (5.1)  is  gs  = 
0.525,  which  is  substantially  smaller  than  the  convergence  rate  per  distributed  relaxation  of  0.9  observed 
with  a  multigrid  FV(2,1)  cycle  applied  to  Eq.  (4.10).  The  smoothing  factor  of  distributed  relaxation  for 
Eq.  (4.10)  is  computed  by  using  local  mode  analysis  as  a  norm  of  the  matrix  G,  defined  as 

G  =  I-(l-gs)(%)-1Mt,  (5.2) 
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Table  2 

The  discretization  errors  in  p  at  convergence  and  the  relative  L\  - norm  errors  after  the  FMG-1  algorithm  for  fully  super¬ 
sonic  flow. 


IN  \‘p 

ns 

HIBm— 

miMiM 

0.8890xl0-4 

1/64 

0.2151xl0-4 

0.02 

1/128 

1/256 

<0.01 

where  matrix  symbols  associated  with  distribution  matrices  for  target  (second-order)  and  driver  (first-order) 
operators  are  denoted  M*  and  Md,  respectively.  The  calculated  smoothing  factor  is  0.752,  appreciably 
greater  than  the  expected  value  of  0.525.  The  increased  smoothing  factor  is  explained  by  the  fact  that 
(M  is  far  different  from  the  identity  matrix. 

With  corrections  distributed  by  the  second-order  distribution  matrix  operator  (Mj  =  M*),  the  smooth¬ 
ing  factor  of  distributed  relaxation  for  Eq.  (4.10)  should  be  identical  to  that  for  relaxation  of  the  Laplacian 
equation.  Implementing  the  defect-correction  relaxation  (5.1)  of  the  Laplacian  equation  but  distributing 
corrections  with  the  second-order  distribution  yielded  an  FV(2,1)  cycle  convergence  rate  per  smoothing  of 
0.51  —  close  to  0.525  predicted  by  the  mode  analysis  of  the  defect-correction  relaxation  (5.1). 

Results  for  an  FMG-1  algorithm,  denoting  a  full  multigrid  algorithm  with  one  FV(2,1)  cycle  at  each 
level,  are  shown  for  fully  subsonic  flow  in  Table  1.  Ail  of  the  norms  used  in  this  and  the  following  tables  are 
Li-norms.  The  discretization  error  in  pressure  is  defined  pointwise  as  e d  =  pj  -p(zj)exact,  where  p(x)exact 
is  the  steady-state  solution  of  Eq.  (3.1)  and  pj  is  the  discrete  solution  of  Eq.  (4.1).  A  second-order  spatial 
convergence  in  the  Li-norm  of  ed  is  evident.  The  total  error  in  pressure  is  defined  as  et  =  p{xj)ex*ct  - pj , 
where  pj  is  the  approximate  solution  obtained  in  the  FMG-1  algorithm.  Small  relative  difference  between 
L\ -norms  of  the  total  and  discretization  errors  (see  Table  1)  indicates  that  optimal  efficiency  has  been 
attained.  The  asymptotic  convergence  rate  of  the  FV(2,1)  cycle  was  0.50-0.52  per  relaxation  for  any  of  the 
four  grids. 

5.2.  Fully  Supersonic  Flow.  In  fully  supersonic  flow,  it  is  generally  possible  to  construct  marching 
methods  to  solve  the  full-potential  equation.  For  first-order  differencing,  the  one-sided  operator  can  easily 
be  solved  over  the  domain.  For  the  second-order  operator  considered  here,  the  defect-correction  relaxation 
with  the  first-order  upwind  driver  can  be  defined  as 

^[{6w)j-2  -  2(6w)j-i  +  (6w)j]  =  -6xx(wj).  (5.3) 

This  scheme  can  be  solved  by  marching  but,  unfortunately,  its  amplification  factor  is  larger  than  1  for  some 
error  components.  A  suitable  full-potential  driver  operator  to  replace  the  left-side  operator  of  Eq.  (5.3)  can 
be  found  in  the  form 


fe+haS^]3,  (5.4) 

where  the  backward  difference  operator  is  defined  to  be  first-order  in  the  above  expression.  This  is  a  5- 
point  operator  that  can  be  solved  by  marching.  The  maximum  amplification  factor  of  the  supersonic  defect- 
correction  relaxation  with  the  driver  (5.4)  and  a  -  0.23  is  0.55.  With  the  second-order  distribution  operator, 


10 


Table  3 

The  discretization  errors  in  p  at  convergence  and  the  relative  L\-norm  errors  after  the  FMG-1  algorithm  for  transonic 
flow  without  a  shock. 


IB 

unSiBl 

0.1075xl0-3 

1/64 

0.2666xl0-4 

1/128 

0.6635xl0-5 

0.03 

1/256 

0.1655xl0-5 

0.11 

Table  4 

The  discretization  errors  in  p  at  convergence  and  the  relative  Li-norm  errors  after  the  FMG-1  algorithm  for  transonic 
flow  with  a  shock. 


h 

INI  :  P 

1M 

:p 

1/32 

0.4396xl0~2 

<0.01 

1/64 

0.2218xl0-2 

0.015 

1/128 

O.llllxlO-2 

0.014 

1/256 

0.5558xl0~3 

<0.01 

computations  for  the  fully  supersonic  flow  corresponding  to  Mach  2  at  inflow  converged  asymptotically  at 
this  rate. 

Results  for  the  FMG-1  algorithm  are  shown  in  Table  2.  The  Lx-norm  of  the  discretization  error  in 
pressure  shows  the  expected  second-order  spatial  convergence.  Optimal  efficiency  is  attained  with  the  FMG- 
1  algorithm.  Here,  the  performance  of  the  single-grid  iteration  at  each  level  is  the  same  as  an  FV(2,1)  cycle, 
because  the  full-potential  factor  is  being  solved  (rather  than  just  smoothed)  by  marching.  The  asymptotic 
convergence  rate  per  relaxation  was  0.53-0.56  for  any  of  the  four  grids. 

5.3.  Transonic  Flows.  Results  for  the  FMG-1  algorithm  are  shown  for  a  smooth  transonic  flow  in 
Table  3.  The  Lx-norm  of  the  discretization  error  in  pressure  shows  the  expected  second-order  spatial  con¬ 
vergence.  The  relative  deviations  of  the  Lx -norm  of  the  total  errors  from  the  Lx-norm  of  the  discretization 
errors  are  increased  from  the  fully  subsonic  or  fully  supersonic  levels  but  still  exhibit  optimal  performance. 
The  asymptotic  convergence  rate  per  relaxation  was  0.52-0.56  for  any  of  the  four  grids. 

Results  for  the  FMG-1  algorithm  are  shown  for  the  transonic  flow  with  a  shock  in  Table  4  .  The  Lx-norm 
of  the  discretization  error  in  pressure  shows  a  first-order  spatial  convergence,  because  there  is  no  limiting  of 
the  interpolation  at  the  shock;  i.e.,  maximum  local  errors  occur  at  the  shock  and  do  not  diminish  with  grid 
refinement.  The  relative  deviations  of  the  Lx-norm  of  the  total  error  with  the  FMG-1  algorithm  from  the 
Li-norm  of  the  discretization  error  again  exhibit  optimal  performance.  The  average  convergence  rate  per 
relaxation  over  10  cycles  varied  from  0.35  on  the  coarsest  mesh  to  0.51  on  the  finest  mesh.  Computations 
for  this  case  with  ENO  differencing  are  presented  in  appendix  C  and  show  a  similar  performance  with  a 
monotone  solution  behavior  in  the  shock  region. 

6.  Concluding  Remarks.  A  general  multigrid  framework  for  obtaining  textbook  efficiency  to  solutions 
of  the  compressible  Euler  and  Navier-Stokes  equations  in  conservation  law  form  has  been  discussed.  The 
general  methodology  relies  on  a  distributed  relaxation  procedure  to  reduce  errors  in  regular  (smoothly 
varying)  flow  regions;  separate  and  distinct  treatments  for  each  of  the  factors  (elliptic  and/or  hyperbolic)  are 
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used  to  attain  optimal  reductions  of  errors.  Near  the  boundaries  and  near  shocks,  additional  local  relaxations 
of  the  conservative  equations  are  necessary.  Example  calculations  are  made  for  the  quasi-one-dimensional 
Euler  equations  for  situations  with  fully  subsonic  and  fully  supersonic  flow,  as  well  as  transonic  flows  with 
and  without  a  shock.  All  of  the  calculations  showed  that  the  FMG-1  algorithm  provides  a  very  accurate 
approximation  of  the  exact  solution. 
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Appendix  A.  Conservative  Fluxes. 

The  fluxes  are  constructed  by  using  the  flux-difference  splitting  of  Roe  [15]  and  are  given  as 

Fj+i  =  l[F(Qfl)  +  F(Qt)  -  |A|(Qh  -  Qt)]  (A.l) 

where  A  =  OF /dQ  is  evaluated  with  a  specific  average  of  the  left  and  right  states,  Qi  and  Qr.  in  order  to 
exactly  satisfy  the  jump  conditions  for  a  shock  [15].  The  left  and  right  states  are  constructed  for  first-order 
accuracy  as 


Q  l  =  Qj , 

Q  R  -  Q.7+1. 

and  for  second-order  accuracy  with  Fromm’s  discretization  as 


(A. 2) 
(A.3) 


Q l  =  Qj  +  j(Qj+i  -  Qi-i)> 

Q  R  ~  Qj+1  —  ^  ( Qj  —  Qi+2  y  ■ 


(A.4) 

(A.5) 


The  dissipation  matrix  |A|  =  T  |A|  T-1  where  T,  T-1  are  diagonalizing  matrices,  the  matrix  of  eigenvalues 


A  = 


Ai  0  0 

0  A2  0 
0  0  A3 


and 


(A-6) 


(Ai ,  A2,  A3)t  =  (u  +  c,  u  -  c,  u)T. 


(A.7) 
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The  tilde  superscript  denotes  that  the  eigenvalues  are  limited  away  from  zero  as, 


(Ad  ii 


l-M  if  |  Aj|  >  e, 

(A2  +  (e)2)/(2e)  otherwise, 


(A.8) 


mainly  to  prevent  expansion  shocks  at  sonic  points  for  the  first-order  discretization  accuracy,  but  also  to 
make  a  smooth  transition  through  the  sonic  point  for  the  full-potential  operator.  The  value  e  is  taken  as 
/3max{|Ai|,  |A2|,  |A3|}.  This  “entropy  fix”  is  not  optimal  in  any  sense,  as  larger  /3  values  for  sonic  points 
are  required  for  the  first  order  scheme  on  coarser  grids.  The  value  of  /J  was  nominally  set  at  0.1  and  was 
arbitrarily  increased  to  0.2  for  grids  of  9  points  or  less  in  the  implicit  first-order  Jacobian  matrices. 


Appendix  B.  Distribution  Matrices. 

The  coefficient  matrices  A ±  are  defined  as  the  positive  and  negative  eigenvalue  contributions  to  A, 
where 


Af  ee  (A,  ±  IAiD/2 


and  A  is  written  in  terms  of  its  eigenvalues,  A(Ai,  A2,  A3),  as 


|(Ai  +  A2)  2Jc(Ai-A2) 

^f(Ai  —  A2)  |(Ai  4-  A2) 

.  *(Ai  -  A2)  ^(A1  +  A2-2A3) 


The  upper  2x2  block  of  L  =  A +8X  +  A  can  be  written  as 


0 

0 

A3 


L=  h 

[  pc(t2)  t\ 

and  the  corresponding  matrix  of  cofactors  is 

M=  fl 

-pc(t2)  ti 

where 

fi  =  l[(At  +  A+)<5-  +  (Ar+A2-)<S+], 
<2  =  i[(A+-A+)<5x-  +  (Ar-A2-)J+]. 

The  determinant  of  L  is 

=  (A+  A  t)S-S~ 

+  (Af  A2  +  Af  Aj)  S~ 5+ 

+  (A|“A2-)<5+<5+ 


(B.l) 


(B.2) 


(B.3) 


(B-4) 


(B.5) 

(B.6) 


(B.7) 


14 


Table  5 

The  discretization  errors  in  p  with  ENO  differencing  at  convergence  and  the  relative  L\-norm  errors  after  the  FMG-1 
algorithm  for  transonic  flow  with  a  shock. 


h 

IMI  -P 

ll£il 

1/32 

0.3597xl0-3 

<0.01 

1/64 

0.9711xl0-4 

0.01 

1/128 

0.5133xl0-4 

0.07 

1/256 

0.2665xl0-5 

0.13 

The  operators  for  second-order  upwind-biased  differencing  are  defined  as 
K  sZ(wi)  =  1^2  (~wi-s  +  2wj-2  +  17,"i-i 

—  36  Wj  +  17wj+i  +  2wj+-i  —  Wj+ 3), 


(B.8) 


Sx  Sx  (wj) 


=  YjTj-2  (™j-4  -  10Wj-3  +  31Wj-2 
-  28wj-i  -  Wj  +  6u)j+i  +  Wj+ 2). 


(B.9) 


Appendix  C.  Tiansonic  Shock  —  ENO  Differencing. 


Fig.  2.  Mach  number  distribution  near  the  shock  computed  with  the  FMG-1  algorithm  compared  to  the  exact  discrete 
solution  (h  —  1/256 ). 

Uniform  application  of  the  upwind-biased  interpolations  for  state  variables  of  appendix  A  (Fromm  s 
scheme)  leads  to  oscillations  at  the  shock.  These  oscillations  can  be  eliminated  by  a  limiting  procedure 
to  reduce  locally  the  order  of  approximation  in  the  region  of  the  shock.  Here,  an  ENO  approach  [16]  is 
used  to  prevent  state-variable  interpolations  from  crossing  the  shock.  At  the  shock  interface,  the  left  (right) 
state  variables,  Ql(Qh),  are  found  by  one-sided  second-order  extrapolation;  at  the  first  interface  upstream 
(downstream)  of  the  shock,  the  state  variables,  Ql(Qh),  are  found  by  second-order  central  averaging. 
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Fig.  2  presents  the  Mach  number  distribution  in  the  region  of  the  shock  and  shows  a  nonoscillatory 
behavior  with  a  one-point  representation  of  the  shock.  The  averaged  convergence  rate  per  relaxation  obtained 
in  the  FV(2,1)  cycle  was  ranged  from  0.36  to  0.44  for  different  target  grids.  The  discretization  errors  and 
relative  deviations  of  the  Li-norm  of  the  total  error  with  the  FMG-1  algorithm  are  shown  in  Table  5; 
comparison  with  Table  4  shows  that  the  ENO  discretization  errors  are  smaller  than  the  discretization  errors 
of  the  unlimited  interpolations,  although  first-order  behavior  is  still  found.  In  both  situations,  for  this 
one-dimensional  case,  second-order  accuracy  is  attained  if  the  la-norm  is  restricted  to  regions  either  far 
upstream  or  far  downstream  of  the  shock.  Although  not  shown,  we  note  that  for  this  case,  uniform  second- 
order  accuracy  can  be  found  with  the  ENO  procedure  above  if  the  entropy  fix  is  dropped  at  the  shock; 
the  shock  jump  is  then  reconstructed  identically,  and  a  zero-point  shock  is  recovered.  The  focus  of  this 
investigation  is  on  convergence,  and  both  Fig.  2  and  Table  5  indicate  optimal  efficiency  has  been  attained 
with  the  FMG-1  algorithm. 


